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A concept for a mathematical model of a gas-electric glassmaking furnace is formulated. A system of equa¬ 
tions for calculating external heat transfer is presented. The computational boundary conditions are examined. 
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Glass is made in regenerative furnaces with transverse 
and horseshoe motion of the fuel-combustion products in the 
melting chamber. The furnaces are divided into gas and 
gas-electric according to the method of heating. Glassmaking 
furnaces can also be classified according to their capacity 
and the chemical composition of the glass. It is obvious that 
the heat and mass transfer processes in each particular fur¬ 
nace develop according to their own definite and characteris¬ 
tic way, which is manifested in the quantitative operational 
parameters of a furnace. On a qualitative level the flow of 
thermophysical processes in furnaces is of a general charac¬ 
ter. For this reason, to construct a mathematical model of a 
gas-electric furnace one must start from the basic principles 
of the general theory of furnaces, where the energetics of the 
technological process is the principal and necessary condi¬ 
tion for process flow. 

The energetic essence of furnace operation unified by the 
concept of the thermal operation of furnaces is the most im¬ 
portant part of the science of the furnace design and calcula¬ 
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tion. The energetic essence of the operation of the entire di¬ 
versity of furnaces by type, dimensions, and purpose can be 
described by a limited number of physical models, by study¬ 
ing which the fundamental principles of furnace design can 
be determined. 

The physical model of a glassmaking furnace can be re¬ 
presented as two zones (Fig. 1): the technological process 
zone (TPZ — the melting tank) and the heat generation zone 
(HGZ — the melting chamber). These zones are divided by 
the molten glass surface F and are separated from the sur¬ 
rounding medium by the refractory masonry M. Since a 
glassmaking furnace is a piece of technological equipment, 
the TPZ is considered to be the primary zone and the HGZ 
the secondary zone, whose function is to create definite en¬ 
ergy conditions in the technological process zone. 

In gas glassmaking furnaces (Fig. la), heat appears in 
the melting tank as a result of heat transfer from the melting 
chamber. The processes ensuring heat flow from the HGZ 
into the TPZ are conventionally called determining and the 
furnaces are said to be heat-exchanger furnaces. The inten¬ 
sity of the technological process in furnaces of this type will 
depend on the efficiency of external heat exchange and the 




Fig. 1 . Diagram of a physical model 
of gas (a), electric ( b ), and gas-elec¬ 
tric (c) glassmaking furnaces. 
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utilization of the heat energy transferred into the melting 
tank. The processes on which the distribution and utilization 
of heat in the technological process zone depend are deter¬ 
mined. 

The determining and determined processes can be of dif¬ 
ferent character. It depends on the conditions of heat flow 
into the tank. For example, in electric glassmaking furnaces 
heat appears directly in the technological process zone 
(Fig. 1 b). For this reason, furnaces with this type of heating 
regime are called heat-generator furnaces. 

A mixed heat-generation regime, where the main source 
of heat delivery to the furnace is the chemical energy of hy¬ 
drocarbon fuel (natural gas, liquid fuel), is realized in gas- 
electric glassmaking furnaces (Fig. lc). The heat introduced 
directly into the molten glass by means of electric power 
(5 - 20% of the heat load of a furnace) is considered to be an 
additional source of energy. It is conventionally called an 
auxiliary electric heater (AEH). 

Thus, the conceptual model of the thermal operation of 
glassmaking furnaces can be represented as a collection of 
heat generation processes and external heat-exchange in the 
melting chamber as well as the hydrodynamics and internal 
heat-exchange in the melting tank. The coupling of the exter¬ 
nal and internal heat and mass transfer provides for the phys¬ 
ical integrity of the thermal operation of a furnace. 

A unique feature of glassmaking furnaces is that the sur¬ 
face of the molten glass can be used to couple the determin¬ 
ing and determined thermal processes. The choice is due to 
the fact that melting of the batch and production of the pri¬ 
mary melt occur on the surface of the tank, most of the heat 
is transferred through this surface from the melting chamber 
into the technological process zone, and finally the tempera¬ 
ture distribution on the surface of the molten glass is a princi¬ 
pal boundary condition for the hydrodynamics of the melting 
tank. In other words, the effect of the determining processes 
on glassmaking is the formation of the temperature field of 
the surface of the molten glass. Such an approach, which 
does not distort the physical essence of the thermal operation 
of the furnace, greatly simplifies the mathematical formaliza¬ 
tion of the thermophysical processes involved in glass¬ 
making, which can be represented by a system of mathemati¬ 
cal models of external and internal heat and mass transfer, 
which are coupled by differential specification of the bound¬ 
ary conditions on the molten-glass surface. 

For the external problem these are boundary conditions 
of the second kind. In this case numerical modeling of the 
mechanisms of external heat exchange with different flame 
organization and variable furnace capacity becomes possible. 
The numerically determined temperature field of the melt 
surface objectively reflects the intercoupled character of the 
heat exchange between the heat generation and technological 
process zones. This gives a basis for using the computed tank 
surface temperatures as boundary conditions of the first kind 
when modeling the internal problem. The correctness of such 
an approach to prescribing a very important condition for 
modeling the hydrodynamics and heat exchange in the melt¬ 


ing tank is due to the thermophysical nature of the formation 
of the convective motion of the molten glass in the melting 
tank of a glassmaking furnace. 

Boundary conditions of the first kind remain admissible 
for the adiabatic surfaces of the melting chamber — the in¬ 
lets of the air channels of the burners and the sections of the 
side walls through which the batch is loaded into the furnace. 
At the same time boundary conditions of the third kind, real¬ 
ized by setting the average heat-transfer coefficient and the 
temperature of the surrounding medium, must be adopted on 
the furnace-barrier surfaces. Such an approach makes it pos¬ 
sible to take account of the design features of the masonry of 
individual sections of the melting part of the furnace as well 
as the effect of the furnace capacity and flame organization 
on the energy parameters of external heat exchange. 

The concept proposed for a mathematical model prede¬ 
termines the functional purpose of the models of external and 
internal heat and mass transfer as well as the requirement for 
the numerical schemes used to solve them. The mathematical 
equations of the external problem and their discrete analogs 
are intended for formalization and calculation of the energy 
aspects of glassmaking. The final equations of the mathema¬ 
tical model of the thermal operation of a furnace — the equa¬ 
tions of hydrodynamics and internal heat exchange — are 
self-consistent. For their solution, the external problem need 
only provide information on the surface temperature field of 
the molten glass. Although these equations are a formaliza¬ 
tion of convective heat-exchange between the melt and melt¬ 
ing-tank surfaces bounding its motion, i.e., they also describe 
the laws of heat transfer, the purpose of this block of the 
mathematical model is purely technological. It is intended 
for determining the conditions for a rational distribution of 
heat in the tank that are required to secure the needed furnace 
capacity and molten-glass quality. 

When using the resolvent zonal method to calculate heat 
exchange in a three-dimensional melting chamber of a fur¬ 
nace the differentiated heat flows along the melting tank, 
which are due to endothermal glassmaking processes and 
heat losses through the tank masonry, are taken into account 
in the source term of the nonlinear equations of heat transfer 
and heat balance for the surface zones of the molten glass. 

A strict mathematical description of the thermal opera¬ 
tion of a glassmaking furnace using an AEH requires that the 
electric field in the melt be calculated accurately. At the pres¬ 
ent time it is impossible to perform such a calculation with¬ 
out having sufficiently reliable data on the electromagnetic 
properties of the molten glass and methods for calculating 
the parameters of the electric field and the electric discharge 
between several electrodes. Taking account of only the inte¬ 
gral characteristics of such a discharge does not solve the 
problem, because the equations of electrodynamics do not fit 
well with the resolvent zonal method of calculating radia¬ 
tive-convective heat exchange. One possible way out this si¬ 
tuation could be to include the AEH power in the source 
terms of the equations for the surface zones of the molten 
glass. 
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Fig. 2. Zones adjoining the zth volume zone. 

To construct a numerical model of external heat ex¬ 
change we assume that the gas motion is stationary, while the 
distribution of the mass velocity in the gas volume is de¬ 
scribed by a known function pv (N), where AUs a point in the 
gas volume. Then, the mass flux G, kg/sec, can be deter¬ 
mined for any section of the surface separating the volume 
zones: 

G = (pv) n F, 

where (pv) n is the average value of the projection of the 
mass velocity of the gas on the normal to the surface F. 

The average value of the heat emission coefficient a, 
W/(m 2 ■ K), calculated using the generalized relation in [1] 
is prescribed for the surfaces zones. 

Let us examine the zonal expression for the convective 
heat flux. Let m volume and n surface zones with the total 
number of zones being l = m + n be separated in the system 
under study. For each zth zone the resulting heat flow is com¬ 
prised of the resulting radiation flow Q - ild , kW, and the con¬ 
vective heat flux Qf on , kW. The convective heat flux of the 

surface zones is understood to be the convective heat emis¬ 
sion flow according to Newton’s law. For the volume zones 
this is the sum of the heat transferred by the moving medium 
and the convective heat-emission flow to adjoining surface 
zones. 

The total (resulting) heat flow Q t is 

a=e, rad +er • 

In the present case, in contrast to the problem of radiative 
heat exchange, it is the quantity Q t that must be determined 
for type-I zones and must be prescribed as condition for 
type-11 zones. For example, in the absence of heat release 
in some volume zone it is necessary to set (Q t = 0 or 
g, rad + g, con = 0 in accordance with the fact that the heat 
emitted by this zone in a stationary regime must be compen¬ 
sated by convective delivery of heat from neighboring zones. 
For zones on the surface of the masonry the quantity Q t is the 


flow of heat losses g, loss . For ideal thennal insulation of the 
masonry Q t = 0 and (?, rad + Qf on = 0, i.e., the heat flow de¬ 
livered to the solid surface by convection must be returned 
into the system as a result of heat exchange by radiation. 

Thus, the state of the system with complex heat ex¬ 
change is characterized by the zonal values of the tempera¬ 
tures Tj and the resulting heat flows Q,. To construct a closed 
system of equations the relations describing radiative heat 
exchange must be supplemented by the dependence of the 
convective heat flow for each zth zone on its temperature and 
the temperatures of the adjoining zones. 

For the prescribed values of the heat emission coeffi¬ 
cients g, con for surface zones is determined from the equation 

Qf on = a ; (Tj — 71 )F t for z = 1, (1) 

where T, K, is the average temperature of the zone and j is 
the number of the volume zone bounding the zth surface 
zone. 

For volume zones the desired relation will have a more 
complicated form, since it must take account of the heat 
emission of the products of combustion to the surface zones 
and the convective transfer of heat between the volume 
zones. 

Let us examine the zth volume zone (i = n + 1, ..., /). We 
shall denote the numbers of the surface and adjoining vo¬ 
lume zones as follows (Fig. 2):yj, ; 2 , ... — surface zones; j [, 
j' 2 , ... — volume zones, of which the products of combustion 
enter the zone under consideration. 

We shall denote by G,-.,- the mass flow transferred to the 
zth zone from the j' zone and we shall take account of the 
fact that an additional mass flow can enter the zth volume 
zone through the adjoining surface zones. 

For a glassmaking furnace the additional mass flow for 
the first volume zone in the path of the gases includes the 
fuel mass flow G s , the air mass flow G, a , and the mass flow 

G. p ' c of the products of combustion drawn to the root of the 

flame in the circulation loops. Then the total mass flow en¬ 
tering (and leaving) the zth volume zone will equal 

£(?,,-=Gf+G a +G i p - c . 

y 

Using the notation introduced above, we shall now ex¬ 
press the individual components of the convective heat flow 
for the zth volume zone: 

- heat flow emitted by the gas to the surface zones 

j 

- heat flow obtained by the zth zone due to the motion of 
the gas 

X c r (Tj. -273 )G j . i + c g (T g - 273) Gf + 
j' 

c a (r a -273)G a + c p c (T p c - 273) G. pc ; 






Numerical Model of External Heat Exchange in a Gas-Electric Glassmaking Furnace 


47 


- heat flow lost by the zth zone due to the motion of the 
gas 


- c i (T. — 273) 


f \ 

VG, +G § +Gf +G?' C 

A—t J t l l l 

V J' 


where c is the specific heat capacity of the medium, 
kJ/(kg • K), and T is the temperature of the products of 
combustion entering from outside the system under study, K. 

Summing these three components we obtain the com¬ 
plete expression for the convective heat flow to the zth vol¬ 
ume zone, which can be written in a compact form as a linear 
combination of temperatures: 

er-IsfcW. ( 2 ) 

k 


where g° is the sum of the terms which do not depend on the 
zonal temperatures and g ki is the coefficient of convective ex¬ 
change, determining the contribution of the Ath zone to the 
convective flux of the zth zone. 

For the volume zones (z = n + 1,...,/) the coefficients g kl 
are determined by the following expressions: 


S°i =c g T g G f + c a F a Gf + c T G. p c + 


+ 273 


Y J (c l -c r )G J , l +(c l -c g )Gf 


(Ci -c a )Gf +{c t -c p c )G 


p.c 


0 for k * j,k ^ /, k + z; 


S ki = 


a jFj for k=j ; 

c r G n for k= f> 



f 

\ 

M 

+ 

VG, +G S +Gf +G pc 

Z— t J 1 i 1 i 

_ j 

K J' 



for k = i, 


where j is the number of the adjoining surface zones and 
j' are the numbers of volume exchange zones from which 
gas enters the zth zone. 

The expression (1) for the convection heat-emission flow 
Q- on for the surface zones (z = 1, ..., n) can be represented in 
form of Eq. (2). Then 


u — A. 

gi 



0 for j, k * 

a i F i for k = /'; 

-a i F i for k = z, 


z; 


where j is the number of the adjoining volume zone. 

Using the general form of the equations of the resolvent 
zonal method for radiation problems and replacing in them 


the resulting radiation flows (7/ ad by the differences 
Qi ~Qr n f° r h = h > —» l we obtain: 

for type-I zones 

Qt = £ ( A krf +g k ,T k ) + gi for i = 1, h ; (3) 

k 

for type-I I zones 

£ Wk + g ki T k ) + g?-Qi= 0 for / = /, (4) 

k 

where A ki is the coefficient of radiative exchange between the 
Ath and zth zones; l = l x + Z 2 is the total number of zones; and, 
/j and I 2 are the numbers of type-I and -II zones. 

To solve the problem the temperatures of type-II zones 
must be determined from the system of equations (4) and the 
values found substituted into the relations (3). In order to 
choose a method for solving the system of equations (4) we 
rewrite it, separating on the left-hand sides of the equations 
the terms of the sum which contain the unknown tempera¬ 
tures of the type-11 zones. For this we denote the sum of the 
first /[ terms in the expressions £ {A ki T k +g ki T k ) in terms 

k 

of the parameters z i : 

z i ='tw?+g»T k ). 

k = 1 

The quantities z i , characterizing for each zth zone the 
part of the resulting heat flow that is due to the interaction of 
this zone with the type-I zones, are determined by the pre¬ 
scribed temperatures of these zones and can be calculated be¬ 
forehand. Substituting these quantities into Eq. (4) we obtain 
a system of equations in which the terms containing the de¬ 
sired temperatures of the type-II zones: 

YtW? +g H T k )+z i +g* i -Q t = 0 

k = l l +l 

for z — /| + 1,..., /. (5) 

In summary, a feature of the description of the compli¬ 
cated heat exchange in the melting chamber of a glassmaking 
furnace is that the system of zonal equations (5) is nonlinear 
with respect to the temperatures of type-II zones. Compared 
with the classical zonal method, the system (5) of nonlinear 
equations obtained above is much easier and more conve¬ 
nient to solve numerically. In the first place, the number of 
equations in the system reduces to l 2 (the number of type-II 
zones), and in the second place the desired temperatures 7j 
(/ = /) + 1,..., 1) enter into the zonal equations explicitly. The 
latter circumstance makes it possible to solve the system (5) 
by Newton’s iteration method, which converges rapidly in a 
wide range. 

The initial data for calculating mass transfer between the 
exchange zones as well as the coefficients of convective ex- 












48 


V. Ya. Dzyuzer et al. 


change are aerodynamic contours of the forward motion of 
the gases flows (sequence of passage through the zones) and 
the degree of fuel burnup [1]. This information makes it pos¬ 
sible, in the first place, to evaluate the power of heat release 
due to gas combustion in the j'th volume zone 

where Qj is the amount of heat released in the ith volume 
zone, kW; B is the gas consumption for glassmaking, 
m 3 /sec; Qf' is the lowest working heat-producing capacity 

of the gas, kJ/m 3 ; and, %, is the degree of fuel burnup in the 
ith volume zone, arb. units. 

Taking into consideration the expression (2) for the con¬ 
vective heat flows Qf on , the equation for calculating the 
mass flux G p c of the products of combustion has the form 


where 


Ski 


0 for k ^ j,k ^ f ,k ;; 

a jFj for k = /; 

CrG} c for k =/; 

X a j Fj + Ci (Gf + Gf + G pc ) for k = i, 


.0 _ S _l /^a + q P-c 

i 

Z C /G p ; c ] + c g G g T & + c a G a r a + T X Cj.G 


Si = 273[(C,. (Gf + Gf + G? c ) - c g G g - c a G a - 


P-C 

p.c j' • 

y 


For the surface zones 

gii = a l F i + K i F=(a l+ K i )F- 


cr° =K F T . 

o i i ± ext i 5 


G p . c =B(l + aL 0 ), 

where L 0 is the stoichiometric air consumption for combus¬ 
tion (consumption factor a = 1.1), m 3 /m 3 . 

Then, for all volume zones in the path of the gases except 
the first one (i = n + 1, ..., /), the total heat delivery due to 
gas motion will be 


Qr =lL a J F J T j - T iT a J F j + 

j j 

Z c j' (Gf, + Gf + Gf' c )71, - 

j' 

273 Z c j’(Gf, +Gf< + G?; c )-c i (Gf +Gf +G ; pc )r i 

j' 

273c ,-(Gf +Gf +Gf- c )=Y J Sj k - gii Ti +g,°, 

k 


where 


gf = 273[(c,Gf -c,Gp + (c,Gf -c Gf,) + 


y " j' 

'-p.c 


+ {c t G?'- c f Gy) ]; 


& 


0 for k ^ /, k f ,k ^ z; 
a-F- for k=j\ 
c r {Gf,+G}+Gj; c ) for k=f- 
Y a j Fj + c, {Gf + Gf + G pc ) for k = z. 


where g.. is the coefficient of convective exchange, which 
takes account of heat exchange between the zth zone and the 
solid surface in contact with it; K t is the coefficient of heat 
transfer in the zth surface zone, W/(m 2 • K); and, 7’ cxt , is the 
temperature of the external surface of the masonry in the zth 
zone, K. 

We shall now make use of the zonal equations (5) to con¬ 
struct a system of working equations of the model for surface 
zones of the molten glass and the lining as well as for the 
volume gas zones. Since the unknowns in the system of non¬ 
linear equations of the resolvent zonal method are the tem¬ 
peratures of the zones, the equations corresponding to the 
type-I zones are eliminated from it. 

For the surface zones of the molten glass with boundary 
conditions of the second kind, the nonlinear equation of heat 
transfer and heat balance assumes the form 

X {A u t£ +g ki T k )-g, T. + A..T* - g gm = 0 

k 

for i=l,..., , (6) 

where A u is the coefficient of radiation exchange between the 
z'th zone and the adjoining zones; n l is the number of type-II 
surface zones (of the molten glass); and, Q gm is a source 
term, kW. 

In turn, 


Qgm -Q 1.1 + 01.2 + 0m.m “ 0e > 


For the first zone in the path of the gases 


e? - I F,Tj - T iH “y F ,- <W. + c,a,T, + 

j j 

f \ 

-r p.C 


V c r G?; c T., -273 c q G q + c a G a + Y c r G?, A 

J J J 8 g a a L-i J j' 


\ 


J 


/ 


c, (Gf + Gf + G pc )7j + 213c t (Gf + Gf + G pc ) - 


Z Ski T k -Su T i +S°i, 


where Q t t and Q t , are the heat consumption for glass¬ 
making and heating the molten glass, respectively, kW; 
0mm represents the heat losses through the masonry of the 
melting tank, kW; and, Q e is the AEFI power, kW. 

A unique relation between the resulting heat flow density 
and the internal temperature 7j can be established for the sur¬ 
face zones of the masonry in the melting chamber. For this 
reason, the boundary conditions of the third kind are set by 
the coefficient of heat transfer F, and the temperature T SUI of 
the surrounding medium. 
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Under these conditions Eq. (5) assumes the form 

k 

K t F t (T i - T sm )+z=Q for z = n, + l, n, (7) 

where z ( . is the heat delivered (consumed) due to a type-I 
zone. 

The heat-release power Q l due to fuel combustion is 
known for the volume gas zones. For this reason, for these 
zones the values of the resulting fluxes Q t = - Q\ are given, 
and the temperatures 7) are determined. Thus, the volume gas 
zones are type-II zones, for which 

X(^ 4 +g ki T k yg.Ti+Aj* +g 1 ° +ef +z,.= o 

k 

for i = n + 1,..., m + n. (8) 

As a result, the numerical model of external heat ex¬ 
change in the melting chamber of the gas-electric furnace is 
formalized by the closed system of zonal equations (6) - (8) 
for the temperatures T h T t . 

When solving a system of nonlinear equations by New¬ 
ton’s method, the emissivity of the gas zone, the matrices of 
the proportionality coefficients, the matrices of the genera¬ 
lized proportionality coefficients, the matrices of the resolv¬ 
ing generalized proportionality coefficients, as well as the 
coefficients of radiation and convection exchange are all cal¬ 


culated at each temperature iteration. The criterion of con¬ 
vergence of the iteration process is that two successive tem¬ 
perature fields in each zone coincide with error not exceed¬ 
ing 10~ 4 . 

The algorithm for calculating the coefficients of radia¬ 
tion exchange is presented in [2], 

The numerical model developed by the present authors 
taken together with the mathematical model of the hydrody¬ 
namics of a melting tank [3] can be used to calculate the pa¬ 
rameters of the thermal operation of high-capacity glass¬ 
making furnaces with gas-electric heating. 
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